function [clusterIdx, H_normalized] = AMKKC(Ks, nCluster, k,lambda)
[nSmp, ~, nKernel] = size(Ks);
%************************************************
% init
%************************************************

w = ones(nKernel, 1) / nKernel;
%*********************************************
% compute Laplacian 
%*********************************************
L= constructAffinity(Ks, k);
[Kwz,KS]=finalKernelGlobalAggregation(Ks,L,w.^2,lambda);
myeps = 1e-5;
iter = 1;
objHistory = [];
isConverged = 0;
maxIter = 100;
while ~isConverged && iter < maxIter
    %*********************************************
    % update H via eigs
    %*********************************************
    H = kernelKmeansViaEigs(Kwz, nCluster);
    obj_H = trace(Kwz*(eye(nSmp)-H*H'));
    objHistory = [objHistory; obj_H];       
    %*********************************************
    % update w via quadprog
    %*********************************************
    z = zeros(nKernel, 1);
    HH = H * H';
    IHH = eye(nSmp) - HH;
    IHH = (IHH + IHH')/2;
    for i1 = 1:nKernel
        z(i1) = sum(sum( KS(:, :, i1) .* IHH));
    end
    opts = [];
    opts.Display = 'off';
    w = quadprog(diag(z), zeros(nKernel, 1), [], [], ones(1, nKernel), 1, zeros(nKernel, 1), ones(nKernel,1), [], opts);
    Kwz = kernelGlobalAggregation(KS, w.^2);
    obj_w =trace(Kwz*(eye(nSmp)-H*H'));
    objHistory = [objHistory; obj_w];
    %fprintf('iter = %d, obj_h = %f, obj_w = %f .\n', iter, obj_H, obj_w);
    if iter > 2 && abs((objHistory(iter-1)-objHistory(iter))/objHistory(iter-1)) < myeps
        isConverged = 1;
    end
    iter = iter + 1;
end
objHistory
H_normalized = H ./ repmat(sqrt(sum(H.^2, 2)), 1,nCluster);
clusterIdx = litekmeans(H_normalized, nCluster, 'MaxIter',100, 'Replicates', 50);
end
function Kw = kernelGlobalAggregation(Ks, w)
[nSmp, ~, nKernel] = size(Ks);
Kw = zeros(nSmp);
for i1 =1 : nKernel
    Kw = Kw + Ks(:, :, i1) * w(i1);
end
end

function [Kwz,KS]=finalKernelGlobalAggregation(Ks,L,w,lambda)
[nSmp, ~, nKernel] = size(Ks);
Kwz = zeros(nSmp);
KS= zeros(size(Ks));
  for i2 =1:nKernel
    Kw = zeros(nSmp);
    K=Ks(:, :, i2);
    N=inv(eye(nSmp)+L(:, :, i2)*K)*L(:, :, i2);
%     for i1=1:nSmp
%        for j1=1:nSmp
%         Kw(i1,j1)=(K(i1,j1)-lambda*K(:,i1)'*N*K(:,j1));
%        end
%     end
    Kw=K-lambda*K'*N*K;
    KS(:,:,i2)=Kw;
    Kwz= Kwz + Kw* w(i2);
  end        
end
function [W] = constructAffinity(Ks, k)
[nSmp, ~, nKernel] = size(Ks);
W= zeros(size(Ks));
for i1 =1:nKernel
    S = zeros(nSmp, nSmp);
    Kss = Ks(:, :, i1)- 10^8 * eye(nSmp);
    [~, indx] = sort( Kss, 2, 'descend');
    indx = indx(:, 1:k);
    lIdx = sub2ind([nSmp, nSmp], repmat([1:nSmp]', k, 1), indx(:));
    S(lIdx) = 1;
    S = max(S, S');
    S= S.* Kss;
    Lw = diag(sum(S,2)) -S;
    Lw = (Lw + Lw')/2;
    W(:,:,i1)=Lw;
end
end
function [H]= kernelKmeansViaEigs(K,cluster_count)
K = (K+K')/2;
opt.disp = 0;
[H,~] = eigs(K,cluster_count,'LA',opt);
end
